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Three-dimensional (3D) g-state Potts models (g=3, 4, and 5) are studied by the tensor product 
variational approach (TPVA), which is a recently developed variational method for 3D classical lat- 
tice models. The variational state is given by a two-dimensional (2D) product of local factors, and is 
improved by way of self-consistent calculations assisted by the corner transfer matrix renormaliza- 
tion group (CTMRG). It should be noted that no a priori condition is imposed for the local factor. 
Transition temperatures and latent heats are calculated from the observations of thermodynamic 
functions in both ordered and disordered phases. 
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I. INTRODUCTION 

The density matrix renormalization group (DMRG), 
which was invented by White in 1992 [|l|, ||], has been 
apphed to a wide class of one-dimensional (_1D) quan- 
tum systems including quantum spin ladders [S . DMRG 
is also efficient for obtaining thermodynamic functions 
of two-dimensional (2D) classical systems ^, |j. Now a 
technical interest in DMRG is to extend its applicability 
to higher dimensional systems |^] . 

It is worth looking at the variational background in 
DMRG in order to obtain a rough image of DMRG in 
higher dimension. In 1995 Ostlund and Rommer showed 
that DMRG assumes so called 'the matrix product wave 
function', and that a very small numbers of parameters 
are sufficient to obtain a good variational energy It 
is a small surprise that such a construction of varia- 
tional state has been known for long years in the field 
of statistical mechanics of 2D classical lattice models. 
In 1945 Kramers and Wannier introduced a very sim- 
ple matrix product as a variational state for the trans- 
fer matrix of the 2D Ising model Later, the idea of 
constructing variational state from local elements was ex- 
tended by Kikuchi 



ter [llUrl, and Vil 



10 (the cluster approximation), Bax- 



ani |l^ (the correlation length equal- 
ity approach). All these approaches calculate the lower 
bounds of the free energies of a 2D system. They use 
a variational state that corresponds to an effective ID 
statistical system with several adjustable parameters. 

Simply increasing the space dimension by one, we can 
extend such variational formula to 3 dimensions. The 
simplest example is the Kramers- Wannier approxima- 



* Electronic address: 
t Electronic address: 



?endiar@savba.sk 



lishinoOphys.sci.kobe-u.ac.jp 



tion applied to the 3D Ising model by Okunishi and 
Nishino , where the 2D Ising model under the external 
magnetic field is treated as variational state, which has 
only two adjustable parameters. The calculated spon- 
taneous magnetization and transition temperature are 
more precise than those obtained from a former attempt 
to extend DMRG to 3D classical systems |15 . A major 
problem in the Kramers- Wannier approximation is that 
one can not always find out a good functional form of 
variational state intuitively, esp. for models other than 
the 3D Ising model. In order to overcome this problem, a 
numerical self-consistent approach has been introduced, 
which we call 'the tensor product variational approach 
(TPVA)' in the following |l|, |l|l. In TPVA the varia- 
tional state is determined automatically, with no refer- 
ence to a priori information on systems. In this paper we 
briefly review the variational principle and the numerical 
algorithm of TPVA, and discuss the applicability of this 
method via trial calculations for (7 = 3,4,5 Potts models. 

In Sec. II we introduce main features of the new algo- 
rithm from the variational point of view. We focus on 
the self-consistent improvement of the variational state. 
A specific way how to apply the variational method to 
the Potts model is presented in Sec. III. We also pro- 
vide the way how to calculate the internal energy and 
the magnetization. The numerical results are presented 
in Sec. IV. In Sec. V we conclude the main results. 



II. VARIATIONAL APPROACH IN TWO 
DIMENSIONS 

For a tutorial purpose we first explain the way how to 
apply TPVA to the square lattice Potts model. (Later in 
the next section we treat the cubic lattice.) 

Let us consider an infinitely long stripe of the width 
2N on the square lattice, which is nothing but the 2N-\eg 



ladder, and consider the g-state Potts model in this finite 
width region. Fig. |^ shows the transfer matrix T[f7|cr] of 
this system when 2N = 6, where 

[cr] = (0-1,(72, ... ,Cr2Ar) and [ct] = (CTi,CT2, • • ■ ,0-2jv) (1) 

represent adjacent rows of g-state spin variables. Here we 
interpret the Potts model as a special case of so called 'the 
interaction round a face' (IRF) model and construct 
T[(7|(t] as a product of plaquette Boltzmann weights 



2N-1 



2N-1 



T[a\a]= H WB{a,a,+,\<J,a,+,) = [] W^b'I^I'^}' 

i=l 1=1 

(2) 

where we have written the nearest neighbor spin pairs 
(Cj (Tj_|_]^) and {a^^^^^l), respectively, as {a} and {a} for 
the book keeping purpose. Following this index rule, the 
local Boltzmann weight is written as follows 
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where we consider the ferromagnetic case (J < 0) 
throughout this paper. 
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FIG. 1: The transfer matrix T[a\a] in Eq. (g). This is the 
case where 2N = 6 and therefore there are five Boltzmann 
weights Wg from i = 1 to i = 5. 

The variational lower bound for the partition function 
per row is the maximum of the Rayleigh ratio 



J2 <S>[d-]T[a\(T]^[a] 



($|r|*) 



(4) 



where ^[a] and ^'[cr] are arbitrary variational states. 
Since the transfer matrix T in Eq. (|^) is symmetric, we 
assume $[(t] = ^'[(t] in the following. 

TPVA consists of local approximations ^ which 
restrict the form of the variational state 5" [a] into a uni- 
form product of local factors 



2N-1 



(5) 



i=l 



where there are only variational parameters. Figure |^ 
graphically represents '^[a] when 2N = 6. A profit of 
writing the variational state in the product form is that 
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FIG. 2: Graphical expression of the variational state ^[<t] in 
Eq. §). 



the norm of the variational state also has the local prod- 
uct structure 

2N-1 



1=1 

2N-1 



= E n {V{a^a,^,))\ (6) 

which is nothing but a partition function of a ID lattice 
model whose local Boltzmann weight is (F((Tj(jj_|_]^)) . In 
the same manner, the numerator of Eq. (Q) is written as 



2N-1 



(*m*>=E n ^^^^i^wi'^kR^^H'^}, (7) 

[a][a] i=l 

which is also a partition function of an effective 2-leg lad- 
der. As we have graphically represented the variational 
state 5'[o'] in Fig. ||, let us also express (^ITI^") graphi- 
cally in Fig. m 
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FIG. 3: Graphical expression of (vl'ITI'I') in Eq. (|). We have 
used black circles for the spins whose configuration sum is 
taken. 

With the use of the variational state thus defined, the 
variational problem in Eq. (jj) is the same as those used 
by Villani Our aim is to obtain the best local 

factor V{(t} numerically. There are several ways to max- 
imize A^jjj. in Eq. (^) , under the condition that the lattice 
size 2N is sufficiently large Keeping the exten- 

sion to three dimensions in our mind, what we consider 
here is to take the variations of X^^^. with respect to each 
local factor 



(5* 



(8) 



When the system size 27V is large enough, it is sufficient 
to consider the variation with respect to the local change 



y{N) _^ y(N) _^ ^y{N) 



(9) 



3 



at the center of the spin row, since we have treated the 
uniform variational state and since the boundary effect 
is neghgible. After a short calculation from the (local) 
extremal condition mn, 



SX 



we obtain an eigenvalue problem 



0, 



(10) 



^ yW{a} = AF^M (11) 



A(^){(7} 



for the local factor V''^\ The new factor A^^^{a} is 
constructed as 



(12) 



whose graphical representation is shown in Fig. 0. The 
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FIG. 4: The factor A'^' in Eq. (|l|) is constructed by joining 
two ^[cr]s and taking spin configuration sum over all spins 
[a] (the black circles) except for the two central ones {ct} = 
((7jv, CTjv+i) (the white circles). 

matrix S'^^ is defined in the same manner 

xw!^\a\cj} , (13) 

where the spin configuration sum is taken over all black 
circles in Fig. ^. 
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FIG. 5: Graphical representation of B*'^' in Eq. 

Since both A^^) and B^^) are constructed from the 
local factor y, the eigenvalue relation Eq. ( pT|) should be 
solved self-consistently. Thus Eq. (11) is a kind of the 



self-consistent equation. A realistic outline how to solve 
the self-consistent equation is as follows. 



• Start the calculation by setting (arbitrary) (f' num- 
bers of initial values for the local factor V{a,a'). 

• Calculate A^^^ and B(^) from Eqs. (|l|) and (|l|), 

respectively, for sufficient large system size 2N . 

• Substitute A^^\ B^^\ and V^^'> to the left hand 
side of Eq. (O) . Obtain the right hand side by 



and normalize it 



V"{a} = 



A^{a} 



V'{a} 
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(14) 



(15) 



• Create a linear combination V^^^ — V + eV" where 
e is a small parameter of the order of 0.1, and re- 
gard it as an improved local factor. After normal- 
izing l^iow go to the second step and repeat the 
calculation till V reaches its (local) fixed point. 

The small parameter e is introduced in order to stabi- 
lize the convergence of the iterative calculation. For sta- 
tistical models that exhibit a phase transition, the self- 
consistent equation has several stable solutions near the 
transition temperature. They correspond to the disor- 
dered state and to each ordered state. In such a case, 
one can 'target' a desired phase just by imposing a very 
small symmetry-breaking field or by setting the initial 
local factor V{a,a') appropriately. 

The main advantage of the above algorithm is that no 
a priori ansatz is necessary for setting up the variational 
parameters. 



III. EXTENSION TO THREE DIMENSIONS 

It is easy to generalize both the variational relation 
(Eq. (Q)) and the construction of the variational state 
in the product form (Eq. (|^)) to 3D models. We can 
increase the space dimension by replacing the row-spin 
[fj] in Eq. (|]) to a 'layer spin' 
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(16) 

Now the system we are considering is an infinitely large 
3D object of size 2N x 27V x oo. As before, we assume 
that the system size 27V is sufficiently large to investigate 
the bulk limit. The 3D generalization of the row-to- row 



4 



transfer matrix in Eq. (g) is a layer-to-layer transfer ma- 
trix. For the 3D g-state Potts model, the layer-to-layer 
transfer matrix is given by 



(27V-1) (2Ar-l) 

rm<^]= n n ^b''^{-i-}^ 

■1=1 j=i 



(17) 



where the IRF-type local Boltzmann weight is written as 
W^'^{a\cj} = W^[ 1(18) 
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We have used the notation i' ~ i + 1 and j' = j + 1, and 
have represented the plaquette spins as {a}. (See Fig. ||.) 



the numerator (4'|T|^') is that of a two-layer 2D lattice 
model 

(2Af-l) (2Af-l) 

E n n V<^^^^{a}wi''^{a\a}V'^^^Ha}. (21) 

[a][a] i=l j=l 

Since the numerator and the denominator are the par- 
tition functions of effective 2D lattice models one can 
calculate both of them using the corner transfer matrix 
renormalization group (CTMRG), which is a variant of 
DMRG applied to 2D lattice models (20). As a byprod- 
uct of CTMRG, the factor A^^'^^a} and the matrix 
_B^^^){a-|cr} can be calculated |2|]. Also, the variational 
free energy per site (F) can be obtained from CTMRG. 
(Numerical details are reported in Ref. |14, p7[.) 

After we obtain the optimized variational factor 
V(^^){a}, the internal energy E and the magnetization 
M can be calculated from A(^^){cr} and B(^^){o-|cr} 
that are created from the optimized variational factor 
YiNN}_ rpj^g internal energy E per site is equivalent to 
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FIG. 6: The IRF type local Boltzmann weight W^^^''{a\a} of 
the q-state Potts models and the variational factors {a} 
and V'-*-''{cr}. The g-state variables a = 0, 1, . . . , g — 1 are 
located at the edges of the cube. We use the notation {a} 
and {a} for the upper and the lower horizontal plaquettes. 

The 2D generalization of the variational state in Eq. (|^) 
can be obtained in the same manner 

(2^-1) {2N-1) 

m = n n ^^''^i^} 

(2^-1) {2N-1) 

= n n V{a,j ai>j aej, a^j') . (19) 
i=i j=i 

There are variational parameters in the local factor 
We assume that the factor is positionally 

independent and the variational state is uniform. The 
local factor at the center of the system is 

The way how to optimize the local factor T/^*-?) , so that 
it maximizes the Rayleigh ratio X^^^ in Eq. (^), is in prin- 
ciple the same as that for 2D systems. The denominator 
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(20) 



is nothing but a partition function of a 2D lattice model 
whose local Boltzmann weight is equal to , and 



E — — J ^{f^N Nt '^N+l n) ~ J ^{'^N '^N N+l) 

- J SifT^^, ^nn) (22) 
and its statistical average is obtained as follows 

E EV{a}B{a\a}V{a} 

{E) = '-^ , (23) 

E V{a}A{a}V{a} 

where we have dropped the superscript (NN) from V, 
A, and B just for simplicity. The magnetization (M) of 
the g-state Potts model can be calculated from the spin 
expectation value 
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(24) 



together with the definition of the order parameter 
q{Sia, 0)) - 1 



(M) 



9-1 



(25) 



IV. NUMERICAL RESULTS 

We calculate the latent heat of the 3D g = 3, 4, and 5 
Potts models, using the internal energy expectation val- 
ues (E) for both ordered and disordered phases. Here- 
after we set /sg = /ig = 1 and only treat the ferromag- 
netic case J = — 1. The convergence control parameter in 
the self-consistent calculation is chosen as e = 0.1. When 
we obtain the variational state for the ordered phase, 
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FIG. 7: The free energy per site {F) of the g = 3 Potts model 
with respect to the inverse temperature K = l/T. 
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FIG. 9: The free energy per site of the g = 5 Potts modeL 
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FIG. 8: The free energy per site (F) of the g = 4 Potts modeL 



we impose a small symmetry breaking field magnetic 
field) to the system during first several iterations, and 
after that we switch it off. For the CTMRG calculations, 
we kept block spin states m up to the value of 20 [ p2[ , 
which is sufficiently large to obtain the thermodynamic 
functions shown bellow. All thermodynamic functions 
converged after 500 iterations at most even in the close 
vicinity of the transition point. 

First we determined the transition temperature from 
the calculated free energy per site (F) with respect to the 
inverse temperature K = l/T. Since the g = 3 5 Potts 
models exhibit the first-order phase transitions, in a close 
vicinity of the transition point Kt there are two minima 
in the free energy F; one corresponds to the disordered 
phase and the other to the ordered one. It is possible to 
detect both of them by way of solving the self-consistent 
equation starting from different initial conditions for lo- 
cal factors. (When the barrier between the minima is 



low, one of the two phases is often accidentally chosen by 
numerical round-off errors. ) Figs. 0, ^ and ^, respectively, 
show the calculated free energy per site {F) for g — 3, 
4, and 5 cases. The black squares and the white circles 
represent {F) for disordered and ordered phases, respec- 
tively, where the point of intersection of these two curves 
results the transition point i^j . The free energy curves 
were drawn by the least-square fitting of plotted data to 



polynomials. The results are, 
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0.6283 for g = 4, and 
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0.5496 for q : 
= 0.672 for q ■ 



For the case q = 3 the most reliable IVIonte Carlo result 
(as far as we know) is ^ = 0.550565 ± 0.000010 [||, 

and thus k''^"^'^ calculated by TPVA is only 0.18% lower 
than K^^. 
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FIG. 10: The energy per site (E) with respect to K for the 
q = 3 Potts model. 

In Figs. |l^, |l^, and |l^, we have plotted the internal 
energies per site (E) as functions oi K = l/T. The latent 
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FIG. 11: The energy per site {E) with respect to K when 
g = 4. 
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FIG. 12; The energy per site {E) with respect to K when 

Q = 5. 



heat is the energy difference 

Q = E+ -E- 



(26) 



between the ordered and disordered phases. As be- 
fore, we have apphed the least-square fittings to inter- 
polate (or extrapolate) the calculated data towards i?+ 

and E~ at the determined transition point 
These energies E'^ and £'~, respectively, are denoted by 
the cross symbols inside the squares and the circles in 
Figs, y, and |l| The resuhs are Q[«=3l = 0.228, 
Q[«=4r~= 0.619, and QI^^^I = 0.693. For the case 
q = 3, QI'^^^I is 41% larger than a Monte Carlo result 

Q 



[9=3] 
MC 



= 0.16160 ±0.00047. 23 



We finally show the calculated spontaneous magneti- 
zation (M) in Figs. |l^, |l^, and |l^. All numerical results 
thus obtained are summarized in Table I. 
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FIG. 13: The magnetization (M) with respect to the inverse 
temperature K for g = 3. 
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FIG. 14: The magnetization (M) with respect to K for q — 4. 



TABLE I: The numerically obtained transition points Kt and 
the latent heats Q by TPVA for the 3D ferromagnetic q—3, 
4, and 5 state Potts models. The values of the m-state block 
spins are given in the second column. 



q 


rn 


Kt 


Q 


3 


20 


0.5496 


0.228 


4 


20 


0.6283 


0.619 


5 


5 


0.672 


0.693 



V. CONCLUSIONS 

Recently proposed self-consistent method for 3D clas- 
sical systems, the TPVA, has been applied to q=3, 4, 
and 5 state Potts models on the simple cubic lattice. 
Thermodynamic functions such as the free energy, the 
internal energy, and the spontaneous magnetizations are 
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FIG. 15: The magnetization (M) with respect to K for g = 5. 



calculated. The numerical algorithm for solving the self- 
consistent equation is stable at any temperature, if the 
convergence control parameter e is chosen to be equal or 
smaller than 0.1. 
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